Renormalization of the Lattice Boltzmann Hierarchy 
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Is it possible to solve Boltzmann-type kinetic equations using only a small number of particles 
velocities? We introduce a novel techniques of solving kinetic equations with (arbitrarily) large 
number of particle velocities using only a lattice Boltzmann method on standard, low-symmetry 
lattices. The renormalized kinetic equation is validated with the exact solution of the planar Couette 
flow at moderate Knudsen numbers. 
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The lattice Boltzmann (LB) method has met with a 
considerable success in a wide range of fluid dynamics 
problems ranging from turbulent to multi-phase flows [H . 
Recently, much of the attention was focused on the use of 
the LB models for simulation of microflows at moderate 
Knudsen numbers, Kn, the ratio of the mean free path 
to a characteristic flow scale 0, 0, 0, 0, 0, 0, 0]- It is 
understood by now that LB models form a well-defined 
hierarchy 0, El, E3- Each level N > 3 of the LB 
hierarchy is characterized by a thoughtfully chosen set of 
discrete velocities whose number scales as N D , where D 
is the spatial dimension. With increasing level N, the 
LB hierarchy constitutes a novel approximation of the 
classical kinetic theory and has to be considered as an 
alternative to more traditional approaches such as higher- 
order hydrodynamics (Burnett or super-Burnett [131 ]) or 
Grad's moment systems [l4|. One salient feature of the 
LB hierarchy, which is crucial for any realistic application 
and which distinguishes it from traditional approaches, 
is that LB is equipped with relevant boundary condi- 
tions derived directly from Maxwell-Boltzmann's theory 
Q. However, proceeding to higher levels N (a must in 
microflow applications) constitutes an increasingly diffi- 
cult computational problem. 

In this Letter, we solve the problem of simulating LB 
models with large velocity sets on small lattices without 
sacrificing any physics or accuracy. The first step in this 
direction is to realize that lower-order models are nothing 
but closures within the higher-order models. This sim- 
ple yet important observation enables us to formulate the 
renormalized LB models of the lower levels in such a way 
that the additional physics of the higher-order models is 
incorporated. In particular, we show with an example of 
exact solution in the stationary Couette flow that the ac- 
curacy of the most commonly used planar D2Q9 LB can 
be enhanced drastically, without introducing additional 
velocities. Thus, we introduce a new way of increasing 
the accuracy of the LB models without significantly in- 
creasing the computational cost. The methodology de- 
veloped herein can be used to renormalize other compu- 
tational kinetic theories. 



We consider the isothermal LB hierarchy of kinetic 
equations 

dtfi + c ia d a fi = Qi(f), (1) 

where /$ are populations of discrete velocities a, i = 
1, . . . , N D , summation convention is assumed, and Q is 
the collision integral satisfying local conservation of den- 
sity and momentum, and vanishing at the equilibrium 
/ eq , where 

fi = W i [P + + l^A \ C ^dP ~ C s d a (3) j ■ (2) 

Here p = Y,iLi fi is the density, j a = J2iLi c iafi is the 
momentum density, c s is the speed of sound, and we shall 
use units in which c s = 1. The weights Wi and the ve- 
locities Cj are so chosen that at each level N the hydro- 
dynamic limit of the kinetic equation ([I]) at low Mach 
numbers is the Navier-Stokes equation. While the hy- 
drodynamic limit of all the models |T]) is the same at 
each level N, their behavior is markedly different when 
exploring the micro-flow domain. Our goal is to mod- 
ify the lowest-level kinetic equations (JTJ) in such a way 
that the non-hydrodynamic features of the higher-level 
models are correctly captured by lower-order models. 

In order to introduce the main ideas, we consider the 
one-dimensional case. For D = 1, the lowest-order (N — 
3) model with three velocities {— v3, 0, v3} (D1Q3) and 
collision integral in the Bhatnagar-Gross-Krook (BGK) 
form, Qi = (f^ q — fi)/r, with a relaxation time r, can 
be written as a moment system for p, j and pressure 

dtp = -d x j, 

d t j = -d x P, , s s 
d t P=-3d xJ --(P-P c «), 

T 

where P cq = p + (j 2 / p) is the equilibrium value of the 
pressure. Note that when writing equation for the pres- 
sure we have used identity for the energy flux, q = 



2 



^n=i c ifi = 3.7: which appears as a consequence of a 
lattice constraint, cf = 3q. The next (JV = 4) member 
of the LB hierarchy is an off-lattice four-velocity model 
based on the roots of the fourth-order Hermite polyno- 
mial {±a, ±6}, where a = \/3 — -\/6 and & = y3 + \/6 
(_D1<34). Assuming a multi-relaxation time collision in- 
tegral, the corresponding moment system reads: 



dtp 
dtj 

d t P 

d t q 



where a = jp- 



-dxj, 
-d x P, 



1 



-d x q--{P~P^), 

T 

-d x {aP + 0p)-hq-q^), 



(4) 



6,0 



-3 are constants 



of the four- velocity set, and q eq = 3j is equilibrium value 
of the energy flux. We have introduced two relaxation 
times, r and 9, in order to distinguish the relaxation of P 
and q. The moment system can be realized, for example, 
as a quasi-equilibrium kinetic equation T^, 1(| with two 
relaxation times. Note that the equations for {p,j, P} 
are not closed within the system (|4|. 

Now it is easy to see that the D1Q3 model §3§ is a clo- 
sure of the D1Q4 moment system (Q}. Indeed, assuming 



9 -C t, and substituting q ps q 



(0) = „eq 



g eq into the equa- 



tion for pressure, one arrives at a closed sub-system for 
{p, j, P} which is identical to ([3]). Note that, from this 
new angle of view, the aforementioned identity for the 
energy flux, q = 3j in ([3]), appears not as a consequence 
of the lattice constraint but rather as an implication of 
the closure. 

Upon realizing this relation between the higher- and 
lower-level "bare" kinetic equations (p} , it is tempting to 
seek improvements for the closure. The simplest way to 
do this is to rescale the time with r, introduce a smallncss 
parameter 77 = 9/t, and compute the hrst correction, so 
that q = q^ + q^\ where 



-TT)0d x p + Ti](3 - a)d x P. 



(5) 



This is certainly in the spirit of the Chapman-Enskog 
method although note that the system to be closed does 
not consist solely of local conservation laws but also in- 
cludes relaxation. The resulting moment system is equiv- 
alent to a renormalized kinetic equation, 



d t f i +c i d x f i -\ i Tr)d*[(a-3)P+0p] 



1 



own, (6) 



with Azp = 1/6 and Ao = —1/3. Thus, we can realize 
the one-step renormalization (OSR) (O (or any other) as 
a kinetic equation for populations on the three-velocity 
lattice supplemented with a source term. The discrete 
velocities of the renormalized kinetic equation (J6j> are on 
the lattice, so that the discretization in time and space 
is straightforward (see, e. g., [13] )■ This simple example 



already demonstrates the advantages of the renormalized 
lattice kinetic equations. 

Now we shall apply the one-step renormalization to the 
particularly important two-dimensional sixteen-velocity 
model (D2Q16, N = 4). The D2Q16 model is a tensor 
product of the two copies of the D1Q4 model considered 
above, and it (or its analogs) has attracted attention re- 
cently 0, 0, Ell as the first LB model which is capable of 
describing correctly the transient Knudsen regime, unlike 
the standard nine-velocity D2Q9 (N — 3) LB model. 

The set of sixteen moments describing the D2Q16 
model is split into the locally conserved (C), slow re- 
laxing (SV) and fast relaxing (Fg) subsystems 

C = {P,j x jy}, (7) 
St {^xx; ^*yyi ^xyi Qxyy-> Qyxx-> "0}") (8) 

F$ {Qxxxi Qyyyi ^x-s ^Pyi ^xj ^yi <i>h (9) 

where ((?) = Ei=i 

Pap = (CaCp), Qaf3~/ = (CuCpCy), 

V = ((c* - l)(cg - 1)), i> a = ((<£ - 3)c x c y ), 

L a = (c a (c 2 x - 3) (c 2 v - 3)), cf> = (c x c y (cl - 3)( C 2 - 3)). 



The closure of the fast subsystem ©, Fg = F cq , where 



^ (0) ={3^,3^,0,0,0,0,0}, 



(10) 



renders the moment sub-system for the nine moments 
C and 5 T equivalent to the moment system of the D2Q9 
model [l[ . Thus, again, the standard LBGK model on the 
nine-velocity lattice appears as a closure of the higher- 
level theory. The one-step renormalization F g is found 
to be (cf. ©): 



(i) 



3rr)(d a p- d a P aa ), 

xx Jy), 

-0W = -3rr]dy(Q XV y - j x ), 

L« = -3rr/d Q V, 
6™ = 0. 



(11) 



With (fTI]) . it is straightforward to write down the renor- 
malized D2Q9 kinetic equation (cf. ©) and to imple- 
ment the space-time discretization. We do not address 
this here. Instead, in order to clearly see the implication 
of the one-step renormalization (|11| , we consider the ex- 
act solution of the renormalized D2Q9 system in the sta- 
tionary Couette flow, where a fluid is enclosed between 
two parallel plates separated by a distance L. The bot- 
tom plate at y = —L/2 moves with the velocity U\ and 
top plate at y — L/2 moves with the velocity t/jj. The 
solution of the renormalized model proceeds essentially 
along the lines of 0]: First, the steady-state OSR D2Q9 
moment system is integrated under the assumption of 
unidirectional flow and no mass flux through the walls. 
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Second, the boundary conditions are applied to compute 
the integration constants of the solution. This step is 
particularly important: The boundary conditions for the 
OSR D2Q9 model are induced by the boundary condi- 
tions of the D2Q16 model. Namely, when the diffusive 
wall boundary conditions Q are applied to the D2Q16 
model, the result is presented in terms of all the moments 
C, S T and F g . We then replace, F e -> F e (0) + F e {1) , and 
obtain the boundary conditions for the OSR D2Q9 sys- 
tem in terms of the C and S T moments. Application of 
the boundary condition completes the solution. Let us 
introduce the mean free path I — ^/3tc s and the Knudsen 
number Kn = l/L. The x-component of the velocity as 
predicted by the OSR D2Q9 model for any r\ = 8/t is: 



u x (y) = Asinh 



( V 



AU + B(j)AU + U, (12) 



where AU = U2 — U1 is the relative velocity of the plates, 
U = (Ui + 17%) /2 is the centerline velocity, and A and B 
are constants which depend only on Kn and 77: 



B = 



A = 



/J.y/rj + 2\/i tanh 



2 v /r;Kn 



(4Kn + A t) v ^+2( A iK n + tanh ^jy^Kn) 
4KnB 

/x 2 07 cosh (2^) + 2^/xsinh (5^) ' 



(13) 



and /i = a + i)~ 3.076. 

It is striking that for 77 = 1 (0 = r), the result (TT2| and 
(TIB")) becomes identical to the one obtained in [?[ for the 
BGK D2Q16 model. We remind (see HQ) that the bare 
D2Q9 model predicts only a linear velocity profile in this 
problem, u x = (l + 2Kn)~ 1 (y/L)A[/ + U, stripped of the 
nonlinear Knudsen layer at the walls. Quite on the con- 
trary, the renormalized D2Q9 model shows clearly the 
Knudsen layer (first term in (jT2J) ) , which is exactly the 
same as in the D2Q16 model itself. The reason for this 
can be traced to the fact that the renormalization re- 
moves the lattice constraint pertinent to the bare D2Q9 
model, namely Q aaa = 3j a . In the present approach, 
this constraint is recognized as a closure relation Fg ' 
which is then corrected by the first term in (TTT|) . Thus, 
the sense of the renormalization is to dress the bare ki- 
netic equations with non-hydrodynamic modes so that 
they reveal correct behavior at non-vanishing Kn. This 
is indeed much in spirit of the renormalization group 
method [20j for spin-lattice models where renormaliza- 
tion improves on the mean-field approximation to dress 
it with correlations. 

In Fig. [U the value of the velocity slip at the wall re- 
sulting from (fT2|) is compared at various Kn with the 
classical data of Willis [2l| for the linearized Boltzmann- 
BGK equation, and with results obtained with the Direct 
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FIG. 1: Slip velocity at the wall as a function of Knudsen num- 
ber. Line: linearized Boltzmann-BGK model 21|; Filled cir- 
cles: DSMC simulation; Open circles: Standard (bare) D2Q9 
model [H, [7}; Open triangles: One-step renormalized D2Q9 
model ([12}, n = 1. 



Simulation Monte Carlo (DSMC) method [13|. The re- 
sult for the bare D2Q9 model is also plotted for compar- 
ison. It is clear that the agreement for the renormalized 
D2Q9 model remains excellent for large values of Kn, 
and the renormalization leads to a drastic improvement 
of the bare D2Q9 model. We conclude this Letter with a 
number of comments using again the simple D1Q4 model 
for the sake of argument: 



(i) The physical meaning of the renormalization is to es- 
tablish an intermediate level between kinetics and hydro- 
dynamics. This intermediate level happens when the dy- 
namics of q becomes slaved by the dynamics of {p, j, P} 
but the dynamics of P is not yet slaved by the dynamics 
of {p, j}. The hydrodynamic limit of model ((H) assumes 
two smallness parameters, e = t/T and fj, = 0/T where 
T is a flow time scale. Instead, we rearrange it in terms 
of two other parameters, e and rj = 9/t = fj,/e. Note that 
r\ need not be small. 

(ii) Although the simple one-step renormalization is quite 
reliable, a rigorous approach to the non-perturbative 
renormalization can be based on the invariance equation 
0, 



A(q) = d t q 



dq 



'V+|!%' + ||^P) =0. (1.4! 



dj 



A stable fixed point of (| 14[) is a fully renormalized q. Ow- 
ing to a specific feature of the LB hierarchy (linearity of 
propagation) , a way to solve Eq. (fT4|) (and similar equa- 
tions in higher dimensions) is the following: (a) Neglect- 
ing the nonlinearity in P cq , we note that the solution q hn 
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of (fl4]) can be found exactly, following the route of exact 
summation of the Chapman-Enskog expansion 24 , [2|| . 
(b) Once the renormalized linear closure g lln is obtained, 
it can be refined to take into account the nonlinearities. 
Substituting q hu into (TT4")) . we compute the defect of in- 
variancc A lm = A(g lm ). With this, a refinement can be 
written, q « q lla + aA lm , where a can be estimated via a 
relaxation method 23]. 

(iii) Importantly, the simple OSR or non-perturbative 
linear renormalization should be sufficient for most of 
the cases of interest in microflow simulations. In fact, 
the nonlinearity of P cq is mainly responsible for the hy- 
drodynamic behavior of the model (advection term in the 
Navier-Stokes equations), whereas the task of renormal- 
ization is to remove lattice constraints and restore such 
features as Knudsen layers, slip velocity etc. With this, 
the renormalized kinetic equations retaining the full P cq 
are still nonlinear, as in the case of Couette flow consid- 
ered above. 

(iv) As a final remark, in the standard kinetic th eory , the 
one-step renormalization was first introduced in [26( as a 
correction to Grad's moment systems, and received some 
attention after the work 27]. However, this approach 
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cannot compete with the LB method both in terms of 
computational efficiency and (more restrictively) because 
of the lack of well-defined boundary conditions. 

In conclusion, the traditional viewpoint on the LB hi- 
erarchy treats each level separately, without any rela- 
tion across the levels. Here, an alternative viewpoint is 
suggested according to which bare kinetic equations of 
the form (fTJ) on the lower and computationally attrac- 
tive levels appear as closures of the higher-level kinetic 
equations. Based on this, we suggested to renormal- 
ize the low-order LB equations in such a way that the 
physics beyond the standard hydrodynamics is correctly 
reported from the higher levels to the lower levels. We 
demonstrated analytically that the renormalized lattice 
Boltzmann model on a standard velocity set reproduces 
the Knudsen layer in the Couette flow which otherwise is 
possible only with the higher-level models. In this sense, 
the renormalized kinetic equations on standard lattices 
are the LB equations, and not the bare ones, written by 
a plain analogy with kinetic theory. We note that the 
renormalization discussed in this Letter concerns prop- 
agation of non-hydrodynamic effects down the LB hier- 
archy and not a renormalization or s ub-g rid modeling of 
Navier-Stokes' turbulence, as in [28|, Eif . We are, how- 
ever, optimistic that the present methods can also be 
useful in the latter problem. 
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